Results (PhD Chapter 2)
Section 2/2
This series of files compile all analyses done during Chapter 2:
- Section 1 regroups maps.
- Section 2 presents regressions and HMSC results.
All analyses have been done with R 3.6.0.
Click on the table of contents in the left margin to assess a specific analysis
Click on a figure to zoom it
To assess Section 1, click here.
To go back to the summary page, click here.
Human activities considered or the analysis:
- city influence: CityInf
- industries influence: InduInf
- dredging dumping sites: DredSit
- commercial ships mooring site: MoorSit
- rainwater sewers: RainSew
- wastewater sewers: WastSew
- city wharves: CityWha
- industries wharves: InduWha
- fisheries (gear used):
- casiers
- bottom-trawling
- dragues
- mesh nets
- palangres
- fisheries (species targeted):
- dogwhelk (Buccinum sp.): Bucc
- common crab (Cancer irroratus): Cirr
- snowcrab (Chinoecetes opilio): Copi
- nordic shrimp (Pandalus borealis): Pbor
- arctic surfclam (Mactromeris polynyma): Mpol
- american lobster (Homarus americanus): Hame
Workspace preparation
# library(G2Sd)
# library(marmap)
library(pander)
library(raster)
library(reshape2)
library(sf)
library(stats)
library(tidyverse)Here, we use data from subtidal ecosystems (see metadata files for more information)
Only stations that have been sampled both for abiotic parameters and benthic species were included.
The script below includes personnal functions, refined data, parameters for each campaign and global means, sd, se.
1. Exploration of relationships between parameters
This section explores relationships between each pair of parameters or AH distances.
Depth
OM
Gravel
Sand
Silt
Clay
CityInf
InduInf
DredSit
MoorSit
RainSew
WastSew
CityWha
InduWha
2. Hierarchical Modelling of Species Communities
This part is under construction.
set.seed(1)
HMSC_X <- apply(all_variables[X], 2, as.numeric)
HMSC_Y <- apply(sp14[,-1], 2, as.numeric)
HMSC_data <- as.HMSCdata(Y = HMSC_Y, X = HMSC_X, Random = Co2014, scaleX = T)
HMSC_prior <- as.HMSCprior(HMSC_data)
HMSC_param <- as.HMSCparam(HMSC_data, HMSC_prior)
HMSC_model <- hmsc(HMSC_data, HMSC_param, HMSC_prior, family = "gaussian", niter = 10000, nburn = 1000, thin = 10)
mixing <- as.mcmc(HMSC_model, parameters = "paramX")
mixing2 <- as.data.frame(mixing)
boxplot(mixing2)
truth <- as.vector(HMSC_param$param$paramX)
average <- apply(HMSC_model$results$estimation$paramX, 1:2, mean)
CI_025 <- apply(HMSC_model$results$estimation$paramX, 1:2, quantile, probs = 0.025)
CI_975 <- apply(HMSC_model$results$estimation$paramX, 1:2, quantile, probs = 0.975)
CI <- cbind(as.vector(CI_025), as.vector(CI_975))
plot(0, 0, xlim = c(1, nrow(CI)), ylim = range(CI, truth), type = "n", xlab = "", ylab = "", main = "paramX")
abline(h = 0, col = "grey")
arrows(x0 = 1:nrow(CI), x1 = 1:nrow(CI), y0 = CI[,1], y1 = CI[,2], code = 3, angle = 90, length = 0.05)
points(1:nrow(CI), average, pch = 15, cex = 1.5)
points(1:nrow(CI), truth, col = "red", pch = 19)
paramXCITable <- cbind(unlist(as.data.frame(average)), unlist(as.data.frame(CI_025)), unlist(as.data.frame(CI_975)))
colnames(paramXCITable) <- c("paramX", "lowerCI", "upperCI")
rownames(paramXCITable) <- paste(rep(colnames(average), each = nrow(average)), "_", rep(rownames(average), ncol(average)), sep = "")
variationPart <- variPart(HMSC_model, c(rep("hab", 7), rep("met", 9)))
barplot(t(variationPart), legend.text = colnames(variationPart), args.legend = list(y = 1.1, x = nrow(variationPart) / 2, xjust = 0.5, horiz = T))
Ymean <- apply(HMSC_model$data$Y, 2, mean)
R2 <- Rsquare(model, averageSp = FALSE)
plot(Ymean, R2, pch=19)
corMat <- corRandomEff(HMSC_model, cor = TRUE)
averageCor <- apply(corMat[, , , 1], 1:2, mean)
corrplot(averageCor, method = "color", col = colorRampPalette(c("blue", "white", "red"))(200))